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2 .  FOREWORD 


In  this  project,  we  have  tackled  a  number  of  matrix 
problems  which  have  wide  applicability.  For  very  large 
sets  of  data,  it  is  necessary  to  perform  numerical 
calculations  so  that  the  investigator  can  have  confidence 
in  the  numberical  results.  Furthermore,  it  is  desirable 
to  perform  the  computation  in  an  efficient  manner,  even 
when  powerful  computers  are  used.  We  have  attempted  to 
devise  methods  which  satisfy  the  properties  of  efficiency 
and  accuracy. 


3 .  NOT  APPLICABLE 


4 .  NOT  APPLICABLE 


.  Statement  of  the  Problem  Studied  and 
b.  Summary  of  the  Most  Important  Results 

Below  are  abstracts  of  the  papers  published  under  this 

grant. 

"Large-Scale  Geodetic  Least-Squares  Adjustment  by 
Dissection  and  Orthogonal  Decomposition,"  by  Gene  H. 

Golub  and  Robert  J.  Plemmons.  Abstract:  Very  large- 
scale  matrix  problems  currently  arise  in  the  context 
of  accurately  computing  the  coordinates  of  points  on 
the  surface  of  the  earth.  Here  geodesists  adjust  the 
approximate  values  of  these  coordinates  by  computing 
least-squares  solutions  to  large  sparse  systems  of 
equations  which  result  from  relating  the  coordinates 
to  certain  observations  such  as  distances  or  angles 
between  points.  The  purpose  of  this  paper  is  to  suggest 
an  alternative  to  the  formation  and  solution  of  the 
normal  equations  for  these  least-squares  adjustment 
problems.  In  particular,  it  is  shown  how  a  block- 
orthogonal  decomposition  method  can  be  used  in  conjunction 
with  a  nested  dissection  scheme  to  produce  an  algorithm 
for  solving  such  problems  which  combines  efficient  data 
management  with  numerical  stability.  The  approach  given 

ff 

here  parallels  somewhat  the  development  of  the  natural 
factor  formulation,  by  Argyris  et  al.,  for  the  use  of 
orthogonal  decomposition  procedures  in  the  finite-element 
analysis  of  structures.  As  an  indication  of  the  magni¬ 
tude  that  these  least-squares  adjustment  problems  can 
sometimes  attain,  the  forthcoming  readjustment  of  the 
North  American  Datum  in  1983  by  the  National  Geodetic 


Survey  is  discussed.  Here  it  becomes  necessary  to 
linearize  and  solve  an  overdetermined  system  of  approxi¬ 
mately  6,000,000  equations  in  400,000  unknowns  -  a  truly 
large  scale  matrix  problem. 

"Truncated  Newton  Methods,"  by  Stephen  G.  Nash.  Abstract: 
The  problem  of  minimizing  a  real-valued  function  F  of  n 
variables  arises  in  many  contexts.  Most  methods  for  solving 
this  problem  have  their  roots  in  Newton's  method,  i.e. 
they  are  based  on  approximating  F  by  a  quadratic  function 
Q.  If  the  number  of  variables  n  is  large,  then  Newton's 
method  can  be  problematic  since  it  requires  the  compu¬ 
tation  and  storage  of  the  Hessian  matrix  of  second 
derivatives.  Use  of  finite-differencing  and  sparse-matrix 
techniques  has  overcome  some  of  these  problems  but  not  all. 

In  this  thesis,  we  examine  a  flexible  class  of 
methods,  called  truncated-Newton  methods.  They  are  based 
on  approximately  minimizing  the  quadratic  function  Q  using 
an  iterative  scheme  such  as  the  linear  conjugate-gradient 
algorithm.  A  truncated-Newton  algorithm  is  made  up  of 
two  sub-algorithms:  an  outer  non-linear  algorithm  con¬ 
trolling  the  entire  minimization,  and  an  inner  linear 
algorithm  for  approximately  minimizing  Q. 

The  most  important  choice  is  the  selection  of  the 

inner  algorithm.  When  the  Hessian  matrix  is  known  to 

be  positive-definite  everywhere,  then  the  basis  linear 
* 

conjugate-gradient  algorithm  can  be  used.  If  not,  Q 
may  not  have  a  minimum.  We  have  used  the  correspondence 
between  the  linear  conjugate-gradient  algorithm  and  the 


Lanczos  algorithm  for  tridiagonalizing  a  symmetric  matrix 
to  develop  methods  for  the  indefinite  case. 

The  performance  of  the  inner  algorithm  can  be  greatly 
improved  through  the  use  of  preconditioning  strategies. 
Preconditionings  can  be  developed  using  either  the  outer 
nonlinear  algorithm  or  using  information  computed  during 
the  inner  algorithm.  A  number  of  diagonal  and  tridiagonal 
preconditioning  strategies  are  derived  here. 

Numerical  tests  show  that  a  carefully  chosen  trun- 
cated-Newton  method  can  perform  significantly  better 
than  the  best  non-linear  conjugate-gradient  algorithms 
available  today.  This  is  important  since  the  two  classes 
of  methods  have  comparable  storage  and  operation  counts, 
and  they  are  the  only  methods  available  for  solving  many 
large-scale  problems. 

"Algorithms  for  computing  the  sample  variance:  analysis 
and  recommendations,"  by  Tony  F.  Chan,  Gene  H.  Golub, 
and  Randall  J.  LeVeque.  Abstract:  The  problem  of  com¬ 
puting  the  variance  of  a  sample  of  N  data  points  (xi)  may 
be  difficult  for  certain  data  sets,  particularly  when  N 
is  large  and  the  variance  is  small.  We  present  a  survey 
of  possible  algorithms  and  their  round-off  error  bounds, 
including  some  new  analysis  for  computations  with  shifted 
data.  Experimental  results  confirm  these  bounds  and 
illustrate  the  dangers  of  some  algorithms.  Specific 
recommendations  are  made  as  to  which  algorithm  should  be 
used  in  various  contexts. 


"Nonorthogonal  analysis  of  variance  using  a  generalized 
conjugate-gradient  algorithm,”  by  Gene  H.  Golub  and 
Stephen  G.  Nash.  Abstract;  A  method  is  developed  that 
computes  an  exact  nonorthogonal  analysis  of  variance 
using  cell  means.  The  method  is  iterative  and  does  not 
require  that  the  non-orthogonal  design  matrix  be  stored 
or  formed.  At  each  stage  in  the  process,  a  balanced 
analysis  of  variance  problem  must  be  solved.  A 
monotonicity  property  for  the  estimates  of  the  regression 
sum  of  squares  is  derived  that  could  be  used  to  minimize 
iteration  in  hypothesis  testing.  An  application  of  the 
algorithm  to  the  solution  of  analysis  of  covariance 
problems  is  also  given. 

"Imposing  curvature  restrictions  on  flexible  functional 
forms,”  by  A.  Ronald  Gallant  and  Gene  H.  Golub.  Abstract 
A  general  computational  method  for  estimating  the  para¬ 
meters  of  a  flexible  functional  form  subject  to 
convexity,  quasi-convexity,  concavity,  or  quasi-concavity 
at  a  point,  at  several  points,  or  over  a  region  are  set 
forth  and  illustrated  with  an  example. 


Gram-Schmidt  Orthogonalization.  Long  experience  of 
using  Gram-Schmidt  has  given  some  insights  into  the  behavior 
of  that  algorithm,  and  a  manuscript  was  submitted  for 
publication.  These  two  problems  are  studied:  first,  what 
happens  when  nonorthogonal  vectors  are  used  in  the  process, 
and  second,  how  orthogonalization  in  elliptic  norms  and  by 
oblique  projectors  should  be  done. 

Gram-Schmidt  orthogonalization  is  a  crucial  part  of 
several  eigenvalue  algorithms.  For  the  Arnoldi  method,  a 
basis  of  the  Krylov  space  is  orthogonalized,  and  in  other 
algorithms  it  is  important  that  the  residual  is  orthogonal 
to  an  appropriately  chosen  subspace.  (Axel  Ruhe) 

Rational  Krylov  Sequences.  This  is  a  starting  point 
for  the  development  of  an  entirely  new  class  of  algorithms 
for  eigenvalue  problems.  The  algorithms  that  are  now  most 
successful  use  vectors  from  a  Krylov  subspace,  which  has 
the  powers  of  the  matrix  as  a  basis  and  thus  yields  poly¬ 
nomials.  Algorithms  are  being  studied  which  are  based  on 
rational  functions  of  the  matrix  (computed  by  shifting  and 
inverting).  These  ideas  have  been  discussed  in  seminars, 
but  neither  experimental  results  nor  written  reports  are 
available  yet. 

It  is  conceivable  that  a  similar  concept  is  usable  for 
computing  matrix  functions,  and  a  possible  application  to 
ODE  solvers  has  been  discussed  with  Prof.  Germund  Dahlquist. 
(Axel  Ruhe) 


Toda  lattice  methods  for  eigenvalue  computation. 

A  method  which  replaces  an  eigenvalue  problem  by  a  set  of 
ordinary  differential  equations,  and  proceeds  by  applying 
an  ODE  solver  to  those,  has  been  proposed  by  Prof.  Deift 
et  al  at  New  York  University.  That  algorithm  has  been 
tested;  also,  a  few  variants  which  involve  shifting  and 
inversion  have  been  developed.  The  experiments  hitherto 
showed  that,  though  innocent-looking,  the  system  needs 
many  steps  of  an  ODE  solver  in  order  to  reach  the  solution. 
An  interesting  fact  is  that  the  number  of  steps  barely  grows 
with  n,  the  order  of  the  matrix,  but  it  is  doubtful  whether 
there  are  any  problems  where  it  really  performs  better  than 
e.g.  QR.  On  the  other  hand,  by  tracking  its  progress 
documented  by  plots,  one  gains  insights  that  are  of  value 
to  understand  the  behavior  of  tridiagonal  eigenproblems , 
and  possibly  get  assistance  in  developing  other  algorithms. 
(Axel  Ruhe) 
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"Nonorthogonal  Analysis  of  Variance  Using  a  Generalized 
Conjugate-Gradient  Algorithm,"  by  Gene  H.  Golub  and  Stephen  G. 
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